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ABSTRACT 

Gravitational lensing is a potentially powerful tool for elucidating the origin of gamma-ray emission 
from distant sources. Cosmic lenses magnify the emission from distance sources and produce time 
delays between mirage images. Gravitationally-induced time delays depend on the position of the 
emitting regions in the source plane. The Eermt/LAT telescope continuously monitors the entire sky 
and detects gamma-ray flares, including those from gravitationally-lensed blazars. Therefore, temporal 
resolution at gamma-ray energies can be used to measure these time delays, which, in turn, can be used 
to resolve the origin of the gamma-ray flares spatially. We provide a guide to the application and Monte 
Carlo simulation of three techniques for analyzing these unresolved light curves: the Autocorrelation 
Function, the Double Power Spectrum, and the Maximum Peak Method. We apply these methods to 
derive time delays from the gamma-ray light curve of the gravitationally-lensed blazar PKS 1830-211. 
The result of temporal analysis combined with the properties of the lens from radio observations yield 
an improvement in spatial resolution at gamma-ray energies by a factor of 10000. We analyze four 
active periods. For two of these periods, the emission is consistent with origination from the core and 
for the other two, the data suggest that the emission region is displaced from the core by more that 
^ 1.5 kpc. For the core emission, the gamma-ray time delays, 23 ± 0.5 days and 19.7 ± 1.2 days, are 
consistent with the radio time delay 26^5 days. 

Subject headings: Galaxies: active - gravitational lensing: strong -gamma-rays: jets -methods: signal 
processing 


1. INTRODUCTION 

Our ability to study gamma-ray radiation from distant 
sources is observationally limited by the poor 0 . 1 °) 
angular resolution of the detectors. We demonstrate a 
technique for using a cosmic lens to enhance the angular 
resolution at gamma-ray energies. 

1 . 1 . History and Challenges 

The high-energy sky is dominated by the most extreme 
and puzzling objects in the universe. These sources 
harbor powerful jets, the largest particle accelerators in 
the universe; they produce radiation ranging from radio 
wavelengths up to very high-energy gamma rays. Obser¬ 
vations of resolved jets reveal very complex structures in¬ 
cluding bright knots, blobs and filaments with sizes rang¬ 
ing from sub-pc up to dozens of kpc (Biretta et al. 1991; 
Siemiginowska et al. 2002; Marshall et al. 2005; Harris & 
Krawczynski 2006; Tavecchio et al. 2007; Marscher et al. 
2008). Observations with improved angular resolution 
reveal that the variability of these sources is very com¬ 
plex; the flaring emission originates from regions close to 
the core, or from knots along the jet, as in M87 (Harris 
et al. 2003, 2006). 

At gamma-ray energies, the Large Area Telescope on¬ 
board the Fermi mission (Eermi/LAT) continuously de¬ 
tects flares from hundreds of blazars. Eermi/LAT scans 
the entire sky in a few hours with excellent time resolu¬ 
tion. However, the angular resolution of Eermi/LAT is 
at best ~ 0 . 1 °, precluding localization of the flares along 
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the jets, even for nearby sources. To resolve gamma-ray 
emission from M87, an improvement in angular resolu¬ 
tion by a factor of 1000 is required. The angular resolu¬ 
tion of pair-creation gamma-ray detectors is typically a 
strong function of gamma-ray energy and viewing angle, 
and is limited by multiple scattering of electron-positron 
pairs and by bremsstrahlung (Atwood et al. 2009). For 
medium-energy gamma-rays, other physical effects, in¬ 
cluding the nuclear recoil, limit the ultimate angular res¬ 
olution of any nuclear pair-production telescope. The 
angular resolution is unlikely to improve by more than 
a factor of a few in future instruments (Bernard et al. 
2014; Wu et al. 2014; Hunter et al. 2014). 

Jets can be well resolved in the radio band down to sub 
parsec scales. The variability time scales at radio wave¬ 
lengths are large, weeks to months, compared to gamma- 
ray time scales of hours to days. Radio monitoring shows 
that roughly 2/3 of the gamma-ray flares coincide with 
the appearance of a new super luminal knot and/or a flare 
in the millimeter-wave core located parsecs from the cen¬ 
tral engine (Marscher et al. 2010; Agudo et al. 2011; 
Marscher et al. 2012; Marscher 2012; Fuhrmann et al. 
2014; Ramakrishnan et al. 2014; Casadio et al. 2015). 

Further study shows correlations and similarities be¬ 
tween multiwavelength and gamma-ray observations 
(Leon-Tavares et al. 2012; Chatterjee et al. 2012; San- 
drinelli et al. 2014). However, Max-Moerbeck et al. 
(2014) modeled the light curves of blazars as red noise 
processes, and found that only 1 of 41 sources with high- 
quality data in both the radio and gamma-ray bands 
shows correlations with a significance larger than 3 (t. 
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They thus demonstrate the difficulties of measuring sta¬ 
tistically robust multiwavelength correlations even when 
the data span many years. 

1 .2. The Importance of the Spatial Origin of Flares 

The origin of gamma-ray flares is a subject of intense 
debate (Nalewajko et al. 2014a; Tavecchio et al. 2010). 
In theoretical modeling, it is generally assumed that 
gamma-ray flares originate from regions close to the cen¬ 
tral engine, typically on parsec scales (Nalewajko et al. 
2014b; Barnacka et al. 2014c; Nalewajko et al. 2014a; 
Hovatta et al. 2015; Tanaka et al. 2011). 

Models where the emission originates from spatially 
distinct knots, as in M87, differ in the underlying funda¬ 
mental physics (Dermer et al. 1992; Bloom & Marscher 
1996; Sikora et al. 1994, 2013, 2009; Giannios et al. 2009; 
Miicke et al. 2003; Bottcher et al. 2013; Reimer et al. 
2004; Stawarz et al. 2005, 2006; Meyer & Georganopou- 
los 2014). Thus, the location of gamma-ray flares is cru¬ 
cial for understanding particle acceleration and magnetic 
fields at both small and large distances from massive 
black holes. 

Multiple variable emitting regions place limitations on 
the use of the most variable quasars for measurement of 
the Hubble parameter based on time delays. Barnacka 
et al. (2015) point out that even a small spatial offset, for 
example 5% of Einstein radius, between the resolved po¬ 
sition of the core and site of variable emission may result 
in a bimodal distribution of values of Hubble parame¬ 
ters characterized by an RMS of ^ 12 [kms“^ Mpc“^]. 
Complex structure can be an important source of sys- 
tematics in measurement of the Hubble parameter from 
gravitationally induced time delays. 

1.3. Gravitational Tensing as a High Resolution Cosmic 

Telescope 

Time delays and magnification ratios derived from 
well-sampled light curves gamma-ray sources can eluci¬ 
date the location of the emitting region. Barnacka et al. 
(2014b) built a toy model placing an M87 analog at 
z ~ 1, with a galaxy acting as a lens at z ^ 0.6. In 
the M87 case, the projected distance between the core 
and the flaring knot HSTl is 60 pc (Biretta et al. 1991) 
or in the toy model ^ 3% of the Einstein radius. In the 
toy model, the differences in time delay (~ 3 days) and 
magnification ratio (^ 0.2) resulting from displacement 
of the radiation source to the jet location can be mea¬ 
sured with current instruments. Thus cosmic lenses are, 
in principle, a tool for revealing the origin of gamma-ray 
flares once the time delays are measured. 

Measurement of the time delays from unresolved light 
curves requires long, evenly-sampled time series with low 
photon noise. Thanks to the Ferm«/LAT observation 
strategy, the light curves are perfectly suited for these 
time delay measurements. Thus, the temporal resolution 
of Fermi/TPTT translates into spatial resolution with the 
help of cosmic lenses. 

1.4. Outline 

We use PKS 1830-211 as a prototypical example of a 
gravitationally-lensed, gamma-ray emitting system (Sec¬ 
tion 2.1;Section 2.2). Our ability to resolve the gamma- 
ray sky relies on time delay measurement (Section 3). 


We evaluate our methods with Monte Carlo simulations 
(Section 3.1). We characterize the signal in Section 3.1.1, 
and evaluate the probability of detecting gravitationally- 
induced time delays and distinguishing them from spuri¬ 
ous fluctuations in the power law noise using approaches 
described in Sections 3.1.2 and 3.1.3. 

We use three methods of time delay estimation: the 
Autocorrelation Function (Section 3.2.1), the Double 
Power Spectrum (Section 3.2.2 and Appendix A), and the 
Maximum Peak Method (Section 3.2.3 and Appendix B). 

We analyze four series of gamma-ray flares in Sec¬ 
tions 4.1, 4.2, 4.3, and 4.4. We use the time delays and 
properties of the lens to elucidate the origin of these flares 
which we discuss in Section 5, and we conclude in Sec¬ 
tion 6. 

2. PKS 1830-211 

PKS 1830-211 is a bright blazar lensed by a galaxy 
located close to the line-of-sight. The emission from 
PKS 1830-211 is detected from radio up to y-ray wave¬ 
lengths. In section 2.1, we describe the gravitationally- 
lensed system and, in section 2.2, we detail the gamma- 
ray emission as observed by the Fermi satellite over a 
6.5-year period. 

2 .1. PKS 1830-211 as Gravitationally-Tensed System 

Pramesh Rao & Subrahmanyan (1988) first recognized 
PKS 1830-211 as a gravitationally-lensed flat spectrum 
radio quasar. The lens is a face-on spiral galaxy at red- 
shift z = 0.886 (Wiklind & Combes 1996, 2001; Winn 
et al. 2002). The quasi-stellar object, a blazar with a 
powerful jet, has redshift z = 2.507 (Lidman et al. 1999). 

The lens distorts the extended radio emission from 
the jet and produces semi-ring-like structure, along with 
two bright images separated by roughly one arcsecond 
(Jauncey et al. 1991). The Australia Telescope Compact 
Array observation at 8.6 GHz revealed these compact 
components, interpreted as emission from the jet core. 

The radio monitoring program lasted for 18 months 
and resulted in the measurement of a time delay of 
26^5 days, and a magnification ratio of 1.52±0.05 (Lovell 
et al. 1998). An independent measurement of the time 
delay using molecular absorption lines yields a delay of 
24^4 days, consistent with the result of radio monitoring 
(Wiklind & Combes 2001). 

The images of PKS 1830-211 have prominent radio 
core-knot structures where the axis is roughly perpendic¬ 
ular to the line of separation of the two mirage images 
of the core. These radio observations define the projec¬ 
tion of the jet. The mirage images of the core are sep¬ 
arated by about 0.98 arcsec and aligned northeast (NE) 
and southwest (SW) in the plane of the sky (Nair et al. 
2005; Garrett et al. 1997; Jones et al. 1996; Guirado et al. 
1999). Jin et al. (2003) reported time dependent posi¬ 
tional variations in the radio centroids at 43 GHz. Marti- 
Vidal et al. (2013) investigated the chromatic variability 
in resolved lensed images of PKS 1830-211 with ALMA 
to probe the jet base of this object. Recently, Marti- 
Vidal et al. (2015) analyzed ALMA data and detected 
a polarization signal (Faraday rotation) related to the 
strong magnetic field at the jet base of PKS 1830-211. 

Winn et al. (2002) modeled the lens parameters. Srid- 
har (2013) refinements yields a best fit for a singular 
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Fig. 1.— Maps based on the lens model of the lensing galaxy for PKS 1830-211 (Section 2.1). Left: Time delays between mirage 
images for different positions of the emitting region within the Einstein radius. The color scale bar represents time delay in days. Right: 
Magnification ratio map for different positions of the emitting regions in the source plane. 


isothermal ellipsoid (SIE) with ellipticity e = 0.091 and 
a lens oriented at 86.1° (for detailed definition of lens 
model see Keeton 2001a). We use these SIE lens pa¬ 
rameters, the redshifts of the source and lens, and the 
Einstein radius of 0.491 arcseconds, to build time delay 
and magnification ratio maps. We base our maps on the 
gravlens code (Keeton 200 la, b). 

Figure 1 (left) shows the time delays between mirage 
images of the emitting region in the source plane. For 
PKS 1830-211, the maximum time delay between the mi¬ 
rage images, which occurs when the source is close to the 
Einstein radius, is ^ 70 days. In the outer region of the 
source plane, the magnification ratio between the mirage 
images > 10. One of the images can be so faint that large 
time delays are undetectable. The right panel of Figure 1 
shows a map of magnihcation ratios between mirage im¬ 
ages of the emitting region located as a function of the 
position in the source plane. 

To determine the position of emitting regions along the 
jet, we need to identify the position of the core and the 
alignment of the jet. Sridhar (2013) determine the posi¬ 
tion of the core based on the time delay (Barnacka et al. 
2011) and magnification ratio reported by Marti'-Vidal 
et al. (2013) along with the Hubble constant (Hinshaw 
et al. 2013). The red circle in Figure 2 delimits the region 
where the core is located. 

Figure 3 shows the predicted total magnification, the 
time delay, and magnification ratios as a function of po¬ 
sition along the jet. Once the time delays and magnifi¬ 
cation ratios are constrained by the light curve, we use 
Figure 3 to identify the positions of the emitting regions. 

2.2. PKS 1830-211 as a Gamma-ray Emitter 
2.2.1. Fermi/LAT Data Analysis 

PKS 1830-211 was detected at gamma rays (Mattox 
et al. 1997; Combi & Romero 1998; Nolan et al. 2012; 
Abdo et al. 2015). The Fermi/LAT telescope is sen¬ 
sitive to photons in the energy range from 20 MeV to 
> 300 GeV (Atwood et al. 2009). Our data analysis 
from 54682 MJD to 57044 MJD, and in the energy range 
from 0.2 — 300 GeV, detects the source at a level of 94 ct. 



x[arcsec] 


Fig. 2.— The range of possible core locations and the jet projec¬ 
tions in the source plane. The gray area shows the allowed range 
(Icr boundary) of the core positions with time delays from 21 to 30 
days (Lovell et al. 1998). The corresponding magnification ratio 
between the resolved images is 1.52 ± 0.05. The blue area rep¬ 
resents the positions of the core constrained by the magnification 
ratio measurement. The red circle delimits the allowed core po¬ 
sitions derived by (Sridhar 2013). Arrows A and B indicate the 
limiting jet projections constrained by resolved radio images. 

We analyze the Fermi/LAT P7REP events and spacecraft 
data using standard likelihood tools distributed with the 
Science Tools v9r32p5 package available at the Fermi 
Science Support Center webpage. 

For source detection and for detection of the time delay, 
we used the P7_CLEAN event selection and the associated 
P7_CLEAN_V6 instrument response function (IRFs). The 
events in the P7_CLEAN class have a high probability of 
being photons. We exclude events with zenith angles 
> 100° from the analysis to limit contamination by Earth 
albedo gamma rays. In addition, we remove events with 
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Fig. 3.— Time delays and magnification ratios as a function of the distance between the emitting region and the core. Left: Total 
magnification defined as the sum of the image magnifications. Middle: Magnification ratios along the limiting jet projections (indicated 
by arrows in Figure 2) Right: Time delays for emitting region located along the limiting jet projections. 



10:00.0 19:00:00.0 50:00.0 40:00.0 30:00.0 20:00.0 10:00.0 18:00:00.017:50:00.0 

Right ascension 


Fig. 4.— Fermi/hAT count map around PKS 1830-211. The 
map contains photons in the energy range from 200 MeV to 300 
GeV. The counts map is smoothed by a Gaussian kernel of cr = 
0.2 deg, with the pixel size of cr = 0.025 deg. 

a rocking angle of > 52° to eliminate time intervals when 
the Earth entered the LAT Field of View (FoV). 

PKS 1830-211 is about ~ 5° from the galactic 
plane (Galactic coordinates, I = 12.°2, b = —5.°7). 
To avoid large background contamination, we ana¬ 
lyze only events with reconstructed energies above 
200 MeV and selected within a square region of 10° 
radius centered on the coordinates of PKS 1830-211 
(i?.A.=278.41333, Zlec=—21.07492). Figure 4 shows the 
count map around PKS 1830-211. The source is clearly 
well-separated from the Galactic plane and there are no 
significant nearby sources. 

To build the gamma-ray light curves, we use a binned- 
maximum likelihood method^ (Mattox et al. 1996). 
This method accounts for sources in the region of in¬ 
terest (ROI), including pulsar PSR J1809-2332. We 
model the background using a galactic diffuse emis¬ 
sion model (gll_iem_v05), and an isotropic component 
(iso_clecin_v05) available at the Fermi Science Support 
Genter webpage. The fluxes are derived from the post- 

^ http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/ 
binned_likelihood_tutorial.html 


launch instrument response functions P7REP_CLEAN_V15. 

The XML source model input to the binned maximum 
analysis contains all the sources included in the Second 
Fermi/LAT catalogue (Nolan et al. 2012) within an an¬ 
nulus of 20° around the ROI. We first analyze the data 
based on the XML source model with free parameters for 
the sources within 7°; the parameters at larger radius 
are fixed to their 2FGL values. We use this XML source 
model to produce the light curves. 

2.2.2. Gamma-Ray Light Curves 

Figure 5 shows the light curve for August 2008 through 
February 2015 with 7 day binning. The energy spectrum 
is well described by a power law with P = 2.54 ± 0.01 
and an integral flux of F(0.2 — 300 GeV) = (1.94 ± 
0.02) X 10“^phcm“^s“^. The highest energy event was 
50 GeV, detected in the time window 55389 MJD - 
55395 MJD. These energies, in principle, are accessi¬ 
ble by the H.E.S.S. II telescope. Thus detection of 
PKS 1830-211 may be possible with H.E.S.S. II. 

Figure 5 shows several active periods. We define active 
periods as times when the gamma-ray emission exceeds 
the average flux by least 2cr. This approach yields four 
active periods. The first series of very bright flares oc¬ 
curs in the period 55420 MJD to 55620 MJD. The sec¬ 
ond series of flares occurs in the period 56050 MJD to 
56200 MJD. Next, a bright single flare occurs around 
July 28, 2014. Recently, on January 8, 2015, another 
flare occurred. Figure 6 shows the light curves of these 
bright flares. 

3. TIME DELAY MEASUREMENT 

Gravitationally-induced time delays are fundamental 
measurements in cosmology. In principle, they provide a 
measurement of the Hubble constant independent of the 
distance ladder (Refsdal 1964; Schechter et al. 1997; Treu 
& Koopmans 2002; Kochanek 2002; Koopmans et al. 
2003; Oguri 2007; Suyu et al. 2013; Sereno & Paraficz 
2014). 

Monitoring of gravitationally lensed sources at both ra¬ 
dio and optical wavelength where the mirage images are 
resolved provide a basis for a number of measured time 
delays (Fassnacht et al. 2002; Eulaers & Magain 2011; 
Rathna Kumar et al. 2013; Tewes et al. 2013; Eulaers 
et al. 2013). Unevenly spaced data resulting from, for 
example, weather and/or observing time allocation, are 
a challenge for light-curve analysis. A number of tech¬ 
niques have been specially developed to utilize these mul¬ 
tiple light curves of mirage images with unevenly sampled 
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Fig. 5.— Fermil'hKJ^ light curve of PKS 1830-211 from August 2008 through February 2015. The fluxes are based on seven-day binning. 
The energy range is 200 MeV to 300 GeV. 



Fig. 6. — Fermi/LAT light curves of flaring PKS 1830-211. We show Flare 1 and Flare 2 with one-day binning (black-filled circles), 
and with 12-hour binning (red open circles). We show Flare 3 and Flare 4 with four-day binning (black-filled circles), and one-day binning 
(red open circles). We show red points only for bins with at least 2cr detection. The green dashed line represents the average flux 
(1.94 it 0.02 X 10“^ photons cm“^ measured from the 7 years light curve of PKS 1830-211 in the energy range from 200 MeV to 300 
GeV. 


data (Edelson & Krolik 1988; Press et al. 1992; Rybicki 
& Press 1992; Burud et al. 2001; Pelt et al. 1998; Pin- 
dor 2005; Scargle 1982; Roberts et al. 1987; Geiger & 
Schneider 1996; Giirkan et al. 2014; Hirv et al. 2011). 

EerTTiz/LAT provides a very long, evenly sampled, light 
curve with low photon noise. The observed light curve 
of lensed blazars is a sum of the mirage images. The 
challenge is to extract the time delay and magnification 
ratio from the time series informed by the model results 
based on shorter wavelength data (Figure 3). 

In the following sections, we investigate three dif¬ 
ferent methods of determining time delays from unre¬ 
solved light curves: the standard Autocorrelation Func¬ 
tion (Section 3.2.1), the Double Power Spectrum method 
(Section 3.2.2), and the Maximum Peak Method (Sec¬ 
tion 3.2.3). Using Monte Carlo simulations, we evalu¬ 
ate the significance levels for these methods, and their 
sensitivity in detecting the gravitationally-induced time 
delays. The Appendices show the detailed steps for the 
Double Power Spectrum (Appendix A.2.1) and for the 
Maximum Peak Method (Appendix B). We use PKS 
1830-211 as a prototype for broader application of these 
techniques. 


3.1. Settings for the Monte Carlo Simulations 

Monte Carlo simulations are a traditional and power¬ 
ful tool for calibrating the analysis of time series. They 
are important in the case of sparsely sampled data and 
they are necessary for evaluating the signihcance of an 
apparent time delay detection (Vaughan 2005). 

In this section, we describe the settings for our Monte 
Carlo simulations. We include the characteristics of the 
time series, the procedures for evaluating the significance 
of time delay detections, and the sensitivity of the meth¬ 
ods for detecting gravitationally-induced time delays in 
unresolved light curves. The full analysis takes advan¬ 
tage of the physical relationship between the time delay 
and the magnification ratio. 

3.1.1. Characteristics of the Signal 

The observed temporal behavior in blazars is rep¬ 
resented by power law noise (Finke & Becker 2014; 
Hayashida et al. 2012; Nakagawa & Mori 2013; 
Sobolewska et al. 2014) where the power spectral den¬ 
sity (PSD) is inversely proportional to the frequency, /, 
of the signal to the power a: 

s{f) cx 1/r ■ 


( 1 ) 
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This random variability is often referred to as noise in¬ 
trinsic to the source (not measurement error), which is a 
result of stochastic processes (Vaughan et al. 2003). As¬ 
tronomers refer to these stochastic fluctuations as signal; 
in other fields, the most common terminology is noise 
(Press 1978). For ease of presentation, we adopt this 
more general terminology and explore properties of var¬ 
ious types of noise. 

Typically, quasars have a ~ 1 — 2 (Equation 1). The 
average slopes for gamma radiation from the brightest 
22 flat spectrum radio quasars and from the 6 bright¬ 
est BL Lacs are 1.5 and 1.7, respectively (Abdo et al. 
2010). During the gamma-ray quiescent state, where 
blazars remain most of the time, the fluctuations in the 
flux are small, and the temporal behavior is characterized 
by power law noise with index ~ 1. During flaring, the 
amplitude of the fluctuation of the flux can increase by a 
few to dozens. In general, the signal is still represented 
by power law noise, but with a greater index a. 

In our simulations, we produce artificial light curves 
with time series represented by red and pink noise. Red 
noise, also known as a Brown noise, has a = 2, consistent 
with the observed behavior of many gamma-ray active 
periods. Figures 6 shows flaring periods of PKS 1830- 
211, and Figure Al shows an example of an artificial light 
curve based on red noise, with and without an artificially 
induced gravitational time delay. The time structure of 
the observed and simulated light curves is remarkably 
similar by eye. 

Pink noise has a power spectrum inversely proportional 
to the frequency of the signal (a = 1); this type of noise 
describes the temporal behavior in the gamma-ray qui¬ 
escent state. For demonstration purposes, we also con¬ 
struct artificial light curves of white noise, a ^ 0. We 
use these types of noise to demonstrate the sensitivity 
of time delay detection to the nature of the underlying 
signal along with the method of analysis. 

We conducted our simulations and analysis using the 
Mat lab environment. We generated samples of power 
law noise using the Little et al. (2007) code. 

The lensed light curve is still power law noise, but it 
contains information about the time delay. The lens it¬ 
self is not a gamma-ray emitter at a detectable level. 
Therefore, we can construct the observed gamma-ray 
light curve as a sum of the lensed components of the 
blazar: 

S(t) = s(t) + s(t + a)/b, (2) 

where S(t) is the unresolved light curve of the lensed 
blazar, composed of the sum of the mirage images. The 
temporal behavior of individual images is determined by 
the source, but the images are shifted in time by the 
gravitationally-induced time delay, a, and with the mag¬ 
nification ratio between mirage images, b. 

We focus on the nature of the gamma-ray emission dur¬ 
ing flaring activity. The durations of these active periods 
range from a few to hundreds of days (see Figure 6). 

The lens model predicts time delays up to ^ 70 days. 
To have a chance of investigating the entire permitted 
range of time delays, the sample has to be at least twice 
as long as the maximum time delay. In our simulations, 
we produced time series of 155 days, exactly the duration 
of the active period of Flare I. 

Fermi/hAT continuously monitors the entire sky, but. 


sometimes, the photon flux of the source can be too low 
to be detected significantly in a time bin chosen a priori. 
There is no gap in the time series, but the flux from 
the source cannot be detected at >2cr. In these cases,we 
compute an upper limit on the flux. The value of upper 
limits does depend, for example, on the exposure within 
the time bin. Thus, if we use upper limits as a measure 
of a flux, satellite-related periodicities can appear in the 
signal. Another approach is to set the flux to zero in 
these low photon flux bins. We check the impact of both 
approaches (Section 4.1). 

One can also interpolate the flux or set the flux to an 
average value. However, this approach could introduce 
a flux, in a particular, that exceeds the observed upper 
limit, thus invalidating the time series analysis. We thus 
do not investigate this approach. 

Monte Carlo simulations provide a strong test of the 
method for treating data points with upper limits. We 
can readily introduce the estimated upper limits into sim¬ 
ulated time series and treat the detections and upper lim¬ 
its in an internally consistent way. This approach gives 
us a measure of the impact of bins with upper limits only 
on the significance of the detection. We discuss this is¬ 
sue further in the results section, where we describe the 
temporal analysis of the PKS 1830-211 light curves. 

3.1.2. Statistical Significance 

A crucial element of the analysis is evaluation of the 
statistical significance of the detection of a time delay. 
An observed modulation of the signal could be a real 
time delay or it could arise purely by chance. Monte 
Carlo simulations provide a way to compute the proba¬ 
bility of detections as opposed to false positives (see e.g. 
Vaughan 2005). 

We produce N = 10® artificial light curves of power 
law noise characterized by a = {0,1,2}. We calculate 
the chance that a particular time delay signal will ap¬ 
pear randomly in the simulated light curve which con¬ 
tains no intrinsic time delays. For each simulated light 
curve, the output consists of time delays and the relative 
power of the signal as a function of the time delay. We 
construct the cumulative probability distribution of false 
detection (CPD: p(POWER)) for each time delay as a 
function of the power. The cumulative probability of a 
false detection is less than p. We investigate p = 0.317 
(1 cr), p = 0.0455 {2 a), p^ 0.0027 (3cr), and p = 0.00006 
(4cr), respectively, for each analysis method and for each 
of the three representative noise spectra. We show the 
results in Sections 3.2.1 and 3.2.2. The CPDs are sen¬ 
sitive both to the underlying noise spectrum and to the 
analysis method. 

3.1.3. Detectability of the Time Delay 

We also use the Monte Carlo simulations to assess the 
chance of detecting a real signal at a given significance 
level. Here the simulated light curves contain artificial 
time delays with the appropriate magnification ratio. Eor 
each combination of noise spectrum, time delay and mag¬ 
nification ratio, the detectability depends on the analysis 
method. In some cases where, for example, the magni¬ 
fication range is large and the time series is too short, 
some methods can not detect the time delay at all. 

The lens model predicts a range of time delays and cor¬ 
responding magnification ratios. For each time delay in 


Resolving the High Energy Universe with Strong Gravitational Leasing 


7 


the range from 1 to 80 days along with the correspond¬ 
ing magnihcation ratios, we produce 10® artificial light 
curves of power law noise with artihcially induced time 
delays. We distribute the time delays uniformly on the 
interval 1 to 80 days with a 1 day interval (the same as 
the light curve binning). These light curves simulate the 
observations of PKS 1830-211. 

We simulate analogous artihcial light curves with red 
noise {a = 2) to evaluate the probability of time delay de¬ 
tection in the observed flaring light curve. To assess time 
delay detectability in quiescent light curves we analyze 
simulated light curves based on pink noise {a = 1). For 
demonstration purposes, we also check the detectability 
of a time delay for white noise light curves (a = 0). 

Analysis of the artihcial light curves gives the proba¬ 
bility of detecting the time delay at 1 cr ,2 ct, 3 a, and 4 a 
levels. From the simulations with no inherent time de¬ 
lays, we know the power where the detection rate of false 
positives is at a given signihcance level. At each power 
(corresponding to a particular signihcance level), we then 
check how many simulated light curves have power above 
this limit. We do this analysis for each time delay and 
signihcance level for false positive detection. 

We discuss the results of this analysis in Sections 3.2.1 
and 3.2.2. As in the case of the analysis of false positives, 
the results for detection of a true time delay depend on 
the input spectrum, the analysis method, and the com¬ 
bination of time delay and magnihcation ratio. 

3.2. Methods 

3.2.1. Autocorrelation Function 

The Autocorrelation Function (ACF) is a standard tool 
for estimating time lags in light curves. The ACF of 
a signal describes the correlation of the values of the 
samples at one time on the values at another. 

For larger a (Eq. 1), the time series contains more ap¬ 
parent structure; there are more and larger amplitude 
variations on shorter time scales. Thus a time delay de¬ 
tection is, in principle, easier for a source characterized 
by larger a. However, the steep spectrum of power law 
noise can also more easily produce spurious peaks in the 
ACF, thus raising the conhdence levels artihcially. 

Figures 7 and 8 show ACF applied to Monte Carlo 
simulations of 10® time series for white, pink, and red 
noise. White noise does not cause spurious time delays 
with large power: the conhdence boundaries are low (^ 
0.3 for Icr), and hat over the full range of time delays 
we explore. However, the probability of detecting time 
delays in the signal characterized by white noise using the 
ACF is very low. The probability of detecting a signal 
at the 1(7 level does not exceed 35%, or 5% at the 2 a 
conhdence level. 

Red and pink noise have a larger potential for produc¬ 
ing spurious peaks in the ACF, (hgure 7); in these cases, 
the conhdence boundaries are elevated. For simulated 
time series, the ACF is most sensitive for time delays in 
the range ^ 5 days to ~ 40 days. The detection of long 
time delays is limited by the hnite length of the sam¬ 
ple; 155 days. Furthermore, mirage images separated by 
longer time delays also have larger magnihcation ratio 
resulting in an attenuation ehect that hinders the detec¬ 
tion. 


3.2.2. Double Power Spectrum 

The Double Power Spectrum (DPS) was the basis for 
detection of the hrst gravitationally-induced time delay 
at gamma-rays (Barnacka et al. 2011). The steps in this 
signal processing are based on widely used methods (Op- 
penheim & Schafer 1975; Brault & White 1971), and 
have been described in Barnacka et al. (2011); Barnacka 
(2013). 

The signal in the time domain is equation (2). The 
Fourier transform of the hrst component, s(t), is s(/), 
and the second component transforms to the frequency 
domain as s(/)e“^’^*'^“. Therefore, the observed signal 
S(t) transforms into 

5(/) = S(/)(l + 6-ie-2"*^“), (3) 

in Fourier space. 

The hrst power spectrum of the source is the square 
modulus of S{f): 

l^(/)P = |S(/)P(1 + b-^ + 2b-^cos{27rfa)). (4) 

The hrst power spectrum is the product of the ’’true” 
power spectrum of the source, |s(/)p, times a periodic 
component with a period (in the frequency domain) equal 
to the inverse of the relative time delay a. Thus, to hnd 
time delay a, we need to hnd the period of the pattern in 
the hrst power spectrum. The Fourier transform of the 
hrst power spectrum, which is in the frequency domain, 
brings us back to the time domain, where the amplitude 
of the signal corresponds to the power of the time delay 
signal present in the original time series. The method is 
similar, in spirit, to the Cepstrum method widely used 
in speech processing and seismology (Bogert et al. 1963). 

We use a Monte Carlo procedure to evaluate the per¬ 
formance of the DPS for power law noise. Figure 9 shows 
the conhdence levels of 1 cr, 2 a, 3 a, and 4 a, for white, 
pink and red noise, respectively. 

The Monte Carlo simulation demonstrates the effec¬ 
tiveness of the signal processing. The DPS method em¬ 
phatically enhances the probability of signihcant detec¬ 
tion of time delays. The probability of detecting a time 
delay 5 — 40 days, at 3 cr level is in the range from 90% to 
40% in contrast with the ACF method where the proba¬ 
bility is 10% at best. The signal processing in DPS elim¬ 
inates the large dependence of signihcance levels and the 
probability of detection of time delay, on the index a of 
power law noise. 

The efficiency of detecting longer time delays is limited 
only by the sample length and the magnihcation ratio be¬ 
tween the mirage images. We follow the same procedure 
independent of the time delay. In principle, one can opti¬ 
mize the signal processing procedure for different ranges 
of time delays to account for the effects of signal atten¬ 
uation and discontinuity at the begin and end of time 
series. 

3.2.3. Maximum Peak Method 

The short duration of the hares relative to the expected 
time delays are an important element in analysis of unre¬ 
solved 7 -ray light curves. Because gamma-ray hares can 
be identihed as distinct events in the time series, and we 
know the range of expected time delays and correspond¬ 
ing magnihcation ratios, we can search for echo hares in 
successive bins directly. 
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Fig. 7. — Confidence level for the autocorrelation method derived from Monte Carlo simulations of 10^ light curves for white (left), 
pink (middle), and red (right) noise. Autocorrelation coefficients as a function of the time delay. The lines indicate Icr, 2(j, 3cr, and 4cr 
confidence levels. The red points show the ACF for Flare 1. 
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Fig. 8 . — Probability of detecting a time delay at a given significance level with the Autocorrelation method. 
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Fig. 9.— Confidence levels for the Double Power Spectrum method. 
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Fig. 10. — Probability of detecting a time delay at a given significance level based on the DPS method. 
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The gamma-ray active state can last for a dozen to 
even hundreds of days; these periods consist of a series 
of individually identifiable flares. Methods like the AGE 
or the DPS are well suited for analyzing these longer 
periods of activity when light curve can be extracted with 
shorter binning ranging from 12 hours to 1 day. When 
the active period consists of a single short flare, bins with 
a longer integration time are necessary. The gamma-ray 
flux before and after these single flares corresponds to the 
quiescent state. Thus to detect the signal at a significant 
level we must increase the size of the bin (exposure). 

The MPM method complements the DPS method for 
these isolated flares. We identify the first brightest flare, 
and calculate the flux ratio between the bin with the 
largest flux (the flare) and flux in successive bins. These 
flux ratios constrain the magnification ratios which are 
not constrained by the DPS method. 

The MPM method enables us to extract additional 
physical constraints from the time series. We compare 
calculated ratios, as a function of a time delay between a 
particular bin and the position of the brightest flare, to 
the magnification ratios as a function of the time delay 
predicted from the model. We identify the time delays 
where the ratio of fluxes is consistent with the predicted 
magnification ratio. These bins might or might not exist. 
If there are bins consistent with the expected delays, the 
data support the picture based on the model. It is impor¬ 
tant to note here that the model takes constraints from 
data at other wavelengths into account. This method 
of analysis of the light curve especially allows us to ex¬ 
clude ranges of time delays where there is no consistent 
magnification ratio observed (~ 80% of the range). 

We demonstrate the method with Monte Carlo simu¬ 
lations in Appendix B. We investigate this approach for 
cases with a single flare and with a series of flares. 

4. RESULTS 

Here, we use PKS 1830-211 as an example of eluci¬ 
dating the spatial origin of the gamma ray flares. In 
particular, we demonstrate that the flares probably do 
not all originate from the same location in the source. 

In the light curve (Figure 5) there are two long active 
periods (red area; Flares 1 and 2) of more than 100 days 
and two isolated individual flares (green area; Flares 3 
and 4). We analyze each of these four flaring periods 
separately. 

4.1. Gamma-Ray Flare 1 

Figure 6 shows Flare 1 with 1-day and 12-hour binning. 
The length of the light curve is 155 days. The temporal 
behavior is characterized by a set of very bright flares. 
Between the flares the flux is close to the long-period 
average (covering the entire light curve); there are upper 
limit detections for ~ 60 days in total. The fit to the 
power spectral density results in a = 1.45, between pink 
and red noise. Figure 11 shows the ACF of Flare 1. 

We investigate two approaches to upper limits; set¬ 
ting an upper limit as the flux, setting the flux to 0 in 
time bins with upper limits. The ACF does not show a 
signihcant difference in time delays and conhdence level 
estimates for these two ways of treating the upper limits. 

The intrinsic variability of the source is consistent with 
the 1 a confidence level. The ACF shows a broad feature 
at a time delay of 17.9 ± 7.1 days at ~ 2cr level. This 



Lags [day] 



Fig. 11.— ACF for Flare 1 along with confidence levels. 

Top: Autocorrelation function for Flare 1 based on upper limits 
as measures of the flux. The confidence levels are based on MC 
simulations of power law noise, with upper limits as measured for 
Flare 1. 

Bottom: Autocorrelation function for the light curve of Flare 1 
with the flux set to 0 in time bins with upper limits. The confidence 
levels are derived by generating time series of power law noise, with 
values set to zero in bins that have measured upper limits. 

result agrees with time delay estimated with the ACF 
performed by Abdo et al. (2015), 19 ± Idays. 

The other broad feature appears at 76 ± 20 days and 
exceeds the 4cr level. However, given the model of the 
lens, this value reaches the maximum allowed time de¬ 
lay. At time delay of ~ 70 days, the magnification ratio 
between the mirage images is larger than 10, and, as we 
have demonstrated with Monte Carlo simulations (Fig¬ 
ure 8), the probability of detecting such gravitationally- 
induced time delays at the 4 ct level using the ACF is 
close to zero. Thus, this feature is probably not produced 
by a gravitationally-induced time delay, but rather re¬ 
flects the time difference between subsets of flares around 
55485 MJD and 55560 MJD. 

Figure 12 shows the analysis of the same time period 
but with the DPS. The DPS method is much more sen¬ 
sitive to signal detection resulting in sharp peaks around 
the time delays. Introducing the values of upper limits 
as a measure of the flux results in a peak at a time delay 
of ~ 52 ± 1.5 day. This time variation corresponds to the 
precession period of the Fermi spacecraft of 53.4 days^. 

This result demonstrates the sensitivity of the DPS in 
detecting even a faint signal in the time series. 

The DPS method, using a 1-day binned light curve, 
detects two time delays at 11±0.5 days and 23±0.5 days 
above the 2cr level. The significance of the detection is 

^ http://fermi.gsfc.nasa.gov/ssc/data/analysis/LAT_caveats_temporal.html 
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Fig. 12.— DPS for Flare 1 in arbitrary units with confidence 
levels. Top: DPS based on upper limits as a measure of the flux. 
Confidence level are based on MC simulations of power law noise, 
including the upper limits as measured for Flare 1. 

Bottom: DPS with the flux set to zero for time bins with upper 
limits. Confidence levels are based on MC simulations also with 
zeros corresponding to upper limits. 



Fig. 13.— The MPM method applied to Flare 1. Red points 
correspond to successive outbursts detected on 55484 MJD. Blue 
points represent the ratio between the outburst at 55560 MJD and 
successive bins. To estimate the errors in the magnification ratio we 
use the flux error at the maxima. The area between the solid and 
dashed lines represents the allowed range of magnification ratios in 
the parameter space defined by the possible projections of the jet 
in the lens plane (Figure 2). The green area represents the range 
of time delays where the observed magnification ratio is consistent 
with the model predictions. 


consistent with expectations for these time delays (see 
Figure 10). For comparison, the DPS method calculated 
for a 12-hour binned light curve yields consistent results, 
with a time delay of 22.5 ± 0.5 days. The time delay at 
11 ± 0.5 days is also present, but at lower significance. 
The lower significance is direct consequence of smaller 
bin. 

To further investigate whether the time delays that ap¬ 
pear in the DPS method are induced by the gravitational 
leasing of a flaring emission region, we use the MPM 
which combines the observations with the predictions of 
the lens model. Figure 13 shows magnification ratios be¬ 
tween the two successive periods following the two largest 
outburst in the Flare 1. To conclude that the detected 
time delay is indeed induced by the gravitational poten¬ 
tial of the lens, we require that both subsets of flares 
have magnification ratios consistent with the time delay. 
The time delay of 11 ± 0.5 days is inconsistent with the 
model. This time delay may be a harmonic of the 23-day 
delay or it may be a previously undetected instrumen¬ 
tal effect. Figure B3 shows an analysis of a randomly 
selected simulated light curve where a harmonic appears 
at this delay. 

The time delay of ~ 23 days is consistent with the 
magnihcation ratio for both subsets of flares. Thus, this 
time delay of 23 ± 0.5 days is probably gravitationally 
induced, and constrains the spatial origin of the Flare 1. 

4.2. Gamma-Ray Flare 2 

Figure 6 shows the light curve for MJD 56043 — 56194. 
The power spectral density is represented by a power 
law with an index a = 1.3. We use this index in Monte 
Carlo simulations to evaluate the confidence levels for 
signal detection. 

The ACF (Figure 7) shows two features at a signifi¬ 
cance level close to 2; the first occurs at a time delay of 
10.1 ± 2.5 days, and the second at 21.1 ± 2.7 days. 

The DPS method, using 1-day binned light curve, 
shows detection of the same features (see Figure 9). The 
first feature appears as a double peak at 11 and 13 days 
at a significance level greater that 2 ct. As in the case 
of Flare 1, the 11-13 day delay is inconsistent with the 
lens model and may be a harmonic or instrumental effect. 
The other peak at 19.7 ± 1.2 days is detected at a signif¬ 
icance level greater then 3 (t. For comparison, the DPS 
method for a 12-hour binned light curve yields consistent 
results, with a time delay of 19 ± 1.0 days. 

For Flare 2, MPM, (see Figure 14), shows a magnifica¬ 
tion ratio consistent with the model predictions for time 
delays in the range from 20 to 23 days. Thus, the time 
delay of 19.7 ± 1.2 days is consistent with the time delay 
expected for the position of the core, and is probably a 
result of gravitational leasing of the flaring gamma-ray 
region. 

4.3. Gamma-Ray Flare 3 

Figure 6 shows the gamma-ray light curve of flare, 
which occurred on July 28. During the flare, the emis¬ 
sion increased by a factor of 5 relative to the average flux. 
The flux for a period of at least 80 days before and after 
Flare 3 is at or below the average flux. 

Figure 15 shows the results of the MPM. The time 
delay range consistent with the expected magnification 
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Fig. 14.— Time delay estimation for Flare 2. Left: Autocorrelation Function with confidence levels. Middle: Double Power Spectrum 
method with confidence levels. Right: Maximum Peak Method, peak 1 corresponds to the flare at 56072 MJD, and peak 2 is a ratio 
calculated relative to the flare at 56146 MJD. Solid and dotted lines indicate predicted magnification ratios along the jet indicated as the 
arrows A and B in Figure 2. 



Fig. 15.— MPM applied to Flare 3 on July 28. 

ratio appears at 46 - 50 days, and corresponds to an in¬ 
crease in the flux recorded at 56912 MJD. To further in¬ 
vestigate whether this period of activity is indeed an echo 
of the flare which occurred at 56865 MJD, we construct 
a light curve around that period with a time binning of 
1 day. 

The delayed counterparts should have similar time evo¬ 
lution. The red points in Figure 15 show that these two 
episodes do not have identical time evolution. The bin 
around 56865 MJD consists of flux close to the average 
for the source; thus the bin may contain significant con¬ 
tribution from the photons originating from the quiescent 
state. 

Flare 3 must have a time delay equal or larger than 48 
days. The secure detection of such a long time delay is 
beyond the sensitivity of the ACF and the DPS for such 
a short light curve with time bins of 4 days. 

If Flare 3 originated from the core, we expect a time 
delay in the range 20 to 30 days. Even with the short 
light curve, we expect to detect the echo flare; in the 
initial flare the flux increases by a factor of ^ 5 and the 
echo flare should appear with a flux at least twice the 
average. Absence of a detection in this range between 20 
and 30 days makes it clear that Flare 3 does not originate 
from the core region. For a time delay > 50 days, Flare 
3 must originate at a projected distance from the core 
> 1.5 kpc (see Figure 3). 

4.4. Gamma-Ray Flare 4 


The most recent gamma-ray activity of PKS 1830-211 
(Figure 6) consists of two flares. For the temporal analy¬ 
sis, we use a light curve with one-day binning consisting 
of 90 days. The ACF, the DPS, and the MPM do not 
show time delays consistent with origination in the core 
(Figure 16). The DPS method indicates a time delay 
at 11.8 ± 0.8 days with a significance of ~ 2a; however 
this time delay is inconsistent with model based on radio 
observations and is thus probably a false positive. The 
time delay of ~ 11 days accidentally corresponds to the 
time between the two flares. The first flare was brighter 
than the average flux for about 4 days and peaked around 
57032 MJD. The second one lasted for about 9 days and 
appeared 2 days after the first one. These flares have 
very different temporal evolution, thus, are not echoes of 
one another. 

Again, the lack of detection of the time delay in the 
range between 20 and 30 days shows that Flare 4 does 
not originate from the core region. The analysis method 
is sensitive for time delays < 50 days and there are no 
other detections. The data show that the time delay 
must be greater than ^ 50 days and thus the radiation 
must originate from a region located at projected dis¬ 
tance from the core > 1.5 kpc (see Figure 3). 

5. DISCUSSION 

Tensing resolves the gamma-ray emission of PKS 1830- 
211 during its flaring periods and limits the origin to the 
core and to regions displaced by > 1.5 kpc along the jet. 
Flares 1 and 2 originate from a region of ~100 pc around 
the core. At the redshift of z = 2.507, where PKS 1830- 
211 is located, a projected distance of 100 pc corresponds 
to ^ 0.02 arcsecond. Thus, this lens improves the angular 
resolution at gamma-ray ~ 10000 times (Figure 17). 

Resolving the high energy universe using cosmic lenses 
relies on the ability to measure time delays and to model 
the mass distribution of the lens. The localization to 
100 pc in this gravitationally-lensed system corresponds 
to an uncertainty in time delay measurement of 5 days. 
The DPS method is an effective approach for measur¬ 
ing the time delay (Section 3.2.2 and Appendix A). This 
method can extract time delays from gamma-ray light 
curves with an accuracy down to 0.5 days. In principle, 
this accuracy can provide a localization of the source to 
^ 10 pc. 

A limiting factor in any lensing analysis is the precise 
model of the lens and alignment of the jet. We have 
used a very conservative position of the core and the jet 
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Fig. 16. — Time delay estimation for Flare 4. Left: ACF with confidence levels. Middle: DPS method with confidence levels. Right: 
MPM, peak 1 corresponds to the bin centered at 57030 MJD, and the peak 2 is a ratio calculated relative to the time bin centered at 
57038 MJD. Solid and dotted lines indicate the predicted magnification ratios along the jet shown as the arrows A and B in Figure 2. 


alignment. In principle, more detailed analysis of all of 
well-resolved radio images can yield better constraints, 
but this analysis is beyond the scope of this paper. 

For these complex sources, there are always puzzles. 
The lack of detection of time delays following Flare 3 
and Flare 4 could indicate some other physical source 
for increased emission. The first possibility is microlens- 
ing. Flare 3 and Flare 4 are unlikely to be microlens- 
ing events because the typical time scale of a caustic 
crossing microlensing event is of the order of months to 
years (Wambsganss 2001); Flares 3 and 4 have a typical 
duration of days and a time structure characteristic of 
gamma-ray flares. 

The size of the emitting region (the source size effect) 
might impact the magnification ratio for Flares 3 and 4. 
A spatially larger emitting region results in a larger mag¬ 
nification ratio (see Figure 2 in Barnacka et al. 2014b). 
However, the minimum variability time scale of ~ 1 day 
observed in these flares constrains the emitting region to 
< 0.01 pc or < 0.001% of the Einstein radius of the lens. 
In other words, the size of the emitting region is small 
enough to have a negligible effect on the magnification 
ratio. 

The final issue is 7 — 7 absorption. Gamma-ray emis¬ 
sion of lensed blazars passes through the lens where low 
energy photons may absorb the gamma rays of one of the 
images passing through the more luminous region of leas¬ 
ing galaxy. The absorption may affect gamma-ray pho¬ 
tons with energies larger than a few GeV. Fermi/LAT 
detects a majority of photons in the energy range > 
100 MeV. In addition, Barnacka et al. (2014a) show that 
the luminosity of a single galaxy is too low to cause sig¬ 
nificant absorption of the gamma-ray flux. If all four 
active periods originated from the same region, absorp¬ 
tion would affect all of them in the same way. However, 
we detect time delays for half of the flaring periods, sug¬ 
gesting that 7 — 7 absorption is irrelevant. 

We have checked the position of the Sun relative to 
the PKS 1830-211. For flares 2,3, and 4, the Sun was 
located outside the Region of Interest (ROI). For flare 1, 
the smallest separation between position of the Sun and 
PKS 1830-211 was ~ 2.5 deg in the period MJD 55558 - 
55561, which is marginal fraction of the Flare 1. Thus, 
the Sun does not affect our analysis. 

A second gravitationally-lensed source, B2 0218-1-35, 
shows behavior similar to PKS 1830-211. The bright 
flaring periods result in time delays consistent with orig¬ 
ination from the core. The time delay measured from 
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Fig. 17.— Resolved positions of gamma-ray flares. The color 
pallet and contours indicate the time delays in days. The long 
arrows show the boundary of the jet alignment limited by well 
resolved radio observations. Gray circle show the position of the 
core from Sridhar (2013). Red circles are further constraints on the 
position of the core using the lens model and the time delay and 
magnification ratio measurements by Lovell et al. (1998). Ellipses 
elucidate the spatial origin of the flares obtained through time delay 
measurement for Flare 1 and Flare 2; they are consistent with the 
core. The top ellipse shows the spatial origin of the time delay 
measured by Barnacka et al. (2011) using the gamma-ray light 
curve in the quiescent state. The short arrow indicates constraints 
from Flares 3 and 4. The time delays > 50 days imply that the 
emitting region must be located at projected distance of > 1.5 kpc 
from the core. 

flaring period at gamma rays is 11.46 ± 0.16 days (Che¬ 
ung et al. 2014). In the radio, Biggs et al. (1999) mea¬ 
sured the time delay of 10.5 ± 0.2 days and Cohen et al. 
(2000) obtained a time delay of 10.1±0.8days. However, 
the most recent gamma-ray flare of B2 0218-1-35 does not 
show delayed counterparts suggesting that in this source, 
flares also have multiple spatial origins. 

6. CONCLUSIONS 

Strong gravitational leasing is a powerful tool for re¬ 
solving the high energy universe. As a prototypical ex¬ 
ample of the power of leasing combined with long, uni- 
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formly sampled light curves in the gamma-ray regime, we 
investigate the spatial origin of flares from PKS 1830-211 
observed with Eermi/LAT. Despite the poor angular res¬ 
olution of gamma-ray detectors, gamma-ray flares can be 
the basis for spatial resolution of a source thanks to the 
unique observational strategy of Eermi/LAT. 

Analysis of four active periods in PKS 1830-211 shows 
that the gamma-ray radiation during two flaring periods 
originated from a region spatially coincident with the 
radio core. The effective spatial resolution we achieve is 
~ 100 pc. 

Two more recent flares apparently do not originate 
from the core region because the time delay must be > 50 
days. This delay and the lens properties derived from 
observations at lower energies indicate that these flares 
must originate at a distance > 1.5 kpc from the massive 
black hole powering the blazar. Our analysis demon¬ 
strates that variable emission can originate from regions 
essentially coincident with the core, or from regions sub¬ 
stantially displaced along the jet. These flares of short 
duration originating from regions within the jet challenge 
our understanding of particle acceleration along with the 
physical conditions along the jet. 

We lay out and apply three methods of time delay 
estimation from unresolved light curves: (1) the Auto¬ 
correlation Function, (2) the Double Power Spectrum, 
and (3) the Maximum Peak Method. We Monte Carlo 
simulations to investigate the strengths and weaknesses 
of these methods and the probability of detecting the 
gravitationally-induced time delays. We provide details 
of their application in the appendices. 

Monte Carlo simulations demonstrate the power of sig¬ 
nal processing (the DPS) in increasing the probability of 
time delay detection over the more standard autocorre¬ 
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APPENDIX 

A. SIGNAL PROCESSING IN THE DOUBLE POWER SPECTRUM METHOD 

The importance of signal processing has been demonstrated in a wide range of applications including analysis of 
speech, imaging, or video and seismic events. Signal processing has played a crucial rule in the development of these 
fields and it is applied extensively in astronomy for spectroscopy, synthetic imaging, and radio astronomy, among 
others (e.g. XCSAO, Kurtz et al. 1992). 

We describe the signal processing method, the Double Power Spectrum, step-by-step. The steps in this signal 
processing are based on widely used methods (Oppenheim & Schafer 1975; Brault & White 1971). We begin by 
preparing the light curve, an input time series. For demonstration purposes, we produce an artificial light curve of 
red noise as a representation of the flaring state. Figure Al shows (red points) the artificial light curve of red noise 
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of duration 155 units. The green light curve shows the same red noise, but including time delays and corresponding 
magnification ratios that simulate gravitational leasing of the source. We induced a time delay of 20 days and a 
magnification ratio of 1.3, following the procedure in Section 3.1.1. 



Fig. A1. — Artificial light curves of red noise. The red points represent pure noise. The green light curve is the sum of two components: 
the first component is the same as the red light curve, and the second component is the identical light curve shifted by 20 days and with 
an a magnification ratio of 1.3. For visualization purposes, we have added error bars of 20% of the average fiux. 




Fig. A2. — No Signal Processing. Right: The First Power Spectrum. Left: The Second Power Spectrum. The SPS is normalized by 
the maximum values in the spectra. The red points represents the FPS and SPS calculated for the red noise light curve (Figure Al). The 
green points show the results for the light curve with an artificially induced time delay of 20 days. The blue points in the SPS indicate the 
result of the full signal processing procedure applied to green light curve that simulates the impact of gravitational lensing. 


The DPS method consists of three major stages. The first is a preparation of the input time series, the light curve. 
The second stage is a calculation of the First Power Spectrum (FPS), and the last stage is a calculation of the Second 
Power Spectrum (SPS). 

Figure Al shows an example of the light curves, used here as the input. For demonstration purposes, we have used 
these light curves and calculated the first and second power spectra, without applying any signal processing. Figure A2 
shows the result of this analysis with the FPS on the left, and the SPS on the right. Note that the light curve (see 
Figure Al) serves as an input to the FPS, the FPS is used as an input for the DPS. 

The green light curve shows the time series with an artificially induced time delay of 20 days. The DPS for this 
signal does not show a hint of the time delay detection when there is no signal processing applied. For comparison, in 
blue, we show the DPS of the same light curve, but after the full signal processing described below. 

A.l. The First Power Spectrum 
A. 1.1. Step 1: Removing the Mean and Windowing 

The input time series in the First Power Spectrum is the light curve shown in Figure Al. We start the signal 
processing by preparing the input. In the first step, we subtract the mean from the time series. This step eliminates 
the large power in the first bin of the first power spectrum. In the next step, we apply a window function to the 
input. Windowing is induced to balance the sharpness of the peak of a periodic signal with the spectral resolution. If 
a time delay is present in the time domain, it will also manifest its presence in the frequency domain (the FPS) as a 
periodic pattern with a period inversely proportional to the time delay (see equation (4)). Thus, we must preserve the 
maximum resolution of the FPS; for this purpose, we use a rectangular window. 
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A. 1.2. Step 2: Zero Padding 

To avoid the large power at low frequencies caused by discontinuity at the beginning and the end of the time series, 
we apply zero padding to the time series. 




Fig. A3.— Right: The First Power Spectrum. Left: The Second Power Spectrum. The input to the FPS are light curve in Figure Al 
after applying Steps 1,2 and 3 of the FPS signal processing. The input to the DPS is the FPS without any signal processing. The colors 
are the same as in Figure A2. 


A.1.3. Step 3: Doubling the Points 

Next, to avoid aliasing, we double the points. Doubling the points does not introduce additional signal, but, it shifts 
the Nyquist frequency and therefore allows the power of the spectrum to go to zero when the frequency approaches 
the Nyquist frequency. We apply the Fourier transform and calculate the power spectrum of the time series. 

Figure Al shows the results of these steps of signal processing. The time delay in the green light curve appears as 
a broad peak around the true value of simulated time delay. 

A.2. The Second Power Spectrum 

The resulting FPS (see Figure A3, left) serves as an input for the SPS. We again process the input before applying 
the Fourier transform. 


A.2.1. Step 1: Flattening and Mean Extraction 

The observed light curve of blazars can be characterized by power law noise. This type of signal has large power at 
low frequencies and thus the signal is not stationary (there is a trend in the data). To ’’flatten” the signal, we take the 
logarithm of the power spectrum (Bogert et al. 1963). Then we can remove the part of the spectrum at low frequencies 
with large amplitude resulting from power law noise. Next, we remove the average from the series. Figure A4 (left) 
shows the input and and the corresponding power spectrum is shown on the right. These steps successfully reduce the 
high amplitudes at short (and spurious) time delays and sharpen the peak around the true time delay. 




Fig. A4.— Right: The First Power Spectrum. Left: The Second Power Spectrum. The spectra are calculated after applying signal 
processing to the light curve, and then processing the FPS by flattening the spectrum and removing a mean (Section A.2.1). 
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A.2.2. Step 2: Windowing and Zero Padding 

In the SPS, the signal of interest is characterized as a peak around the true value of the time delay. Thus, successful 
processing should sharpen the peak. Here we use the Bingham window^, a combination of a rectangular and a Hanning 
window. 




10 20 30 40 50 60 70 80 

Lags [day] 


Fig. A5.— Right: First Power Spectrum. Left: Second Power Spectrum. Results obtained after applying the full signal processing. 

Finally, we apply zero padding to the signal and calculate the final SPS (see Figure B3). Monte Carlo simulation 
demonstrate the effectiveness of the signal processing we have applied. Figure B3 displays the results for one possible 
realization of the light curve. The DPS method is efficient in detecting time delays independent of the character of 
the signal (whether the light curve is white or red noise) and the method is also very resistant to producing spurious 
detection of time delays even in very structured time series of red noise. 

B. MAXIMUM PEAK METHOD - SIMULATIONS 

The time series of gamma-ray light curves consist of short duration flares with variability time scales on the order 
of hours to days. We can use this characteristic of the time series to identify the most prominent flares in the light 
curves; we calculate the ratio between the flux during the flare and the flux in the successive bins. The resulting flux 
ratio is a proxy for the magnification ratio. 

The mirage image arriving first has a larger magnification than the echo images. Thus, the most luminous flares 
can be associated with the first mirage image. The echo flares should appear with a flux diminished by a factor 
corresponding to the magnification ratio. Thus we can look for time periods (time bins) when the flux ratios (normalized 
to brightest flares) are consistent with the magnification ratios predicted for a particular projection of the jet and lens 
model. 

The active periods may consist of a series of flares (Flare 1 and Flare 2), or they can consist of a single outburst 
(Flare 3 and Flare 4). Analysis of light curves containing single flares is difficult because of the low photon statistics 
before and after the flare. Thus the light curve calculated for short time bins may consist primarily of upper limits. Such 
light curves are useless for the Autocorrelation Function, or even for the Double Power Spectrum methods. However, 
light curves containing single isolated flares can still reveal the echo flares and their time delays and magnification 
ratios. We demonstrate the Maximum Peak Method applied to a single simulated flare in Section B.l. 

Long active periods with good photon statistics allow construction of gamma-ray light curves where only a small 
fraction of the bins have upper limits only. These light curves are perfectly suited for methods like the Autocorrelation 
Function and the Double Power Spectrum. In this case, the Maxim Peak Method is a powerful consistency check if the 
detected time delays are consistent with the parameters predicted by lens model. We demonstrate the performance of 
this method on simulated superimposed series of flares in Section B.2. 

The Maximum Peak Method is complementary to the Autocorrelation Function and the Double Power Spectrum in 
the search for gravitationally-induced time delays in unresolved flaring light curves. 

B.l. Single Flares 

Flare 3 and Flare 4 are single isolated outbursts. We base our simulations on Flare 3. 

Flares are observed when the emission significantly exceeds the level typical of the quiescent state. For single isolated 
flares the emission before and after the flare is consistent with the average flux. Therefore, the average flux originates 
from the site of quiescent emission. The temporal behavior of quiescent emission is accurately represented by pink 
noise. We thus start our Monte Carlo simulations by producing light curve composed of pink noise. The flaring episode 
with a flux increase by a factor of 5 is inconsistent with pink noise. We add a flare to the light curve with a time 
structure and flux per each bin similar to Flare 3. Then we introduce an echo flare with a time delay of 48 days 
and magnification ratio 4.5 (Figure Bl, green points). The choice of these parameters relays on the results of the 
Maximum Peak Method for Flare 3. We also simulate an echo flare with time delay of 23 days and magnification ratio 

^ http://www.vibrationdata.com/tutorials/Bingham_compensation.pdf 
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1.8 (Figure Bl, red points). We simulate the echo flare at 23 days to demonstrate that if the flare would originate 
from the region consistent with the core then the echo flare would be detectable. 


x10® 



Fig. Bl.— The artificial light curve generated as pink noise with a flare like structure with time delay of 48 days. This light curve 
simulates Flare 3. We include an echo flare with a time delay of 48 days and magnification ratio of 4.5 (green points). Red points represent 
the light curve with an echo flare at a time delay of 23 days and magnification ratio of 1.8. 


Figure B3 shows the result of applying the Maximum Peak Method to the simulated light curve shown in Figure Bl. 
The method shows that the ratio we obtain between the flux of the flare peak and the flux in the bin corresponding 
to echo flare agree with model predictions. The method rejects the majority of time delay ranges where there is no 
consistent magnification ratio. 




Fig. B2.— Maximum Peak Method applied to the simulated light curve Figure Bl. The red area indicates the time bin corresponding 
to the simulated time delay. Solid and dashed lines indicate the model predictions for magnification ratio as a function of the time delay 
for boundary alignments of the jet constrained by the radio observations (Figure 2). Left: Results for simulated light curve with induced 
echo flare with time delay of 23 days and magnification ratio 1.7. Right: Simulated light curve with induced time delay of 48 days and 
magnification ratio 4.5 


Note that the maximum distance of the emitting region along the jet which we can constraint depends on the ratio 
between the observed flux of the flare and the flux relative to quiescent state. Flare 3 and the simulated flare exceed 
the average flux by a factor of ~ 5. In this example, if the predicted magnification ratio is larger than ^5, we do not 
expect to be able to detect the echo flare; the flux of the echo flare is then below the average of the quiescent state. 
Thus, we can only test expected magnification ratios for Flare-3 in the range from 1 to 5. 

For Flare-3 the magnification ratio ~ 5 corresponds to a region along the jet located at least 1.5 kpc from the core 
(see Figure 3). Detection of a consistent magnification ratio f ^ 5 still does not provide clear evidence of echo flare 
detection because the data are also consistent with a flare from a region at distances > 1.5 kpc from the core. In other 
words, the observed magnification ratio sets a very interesting limit on the distance between the core and the origin 
of the flare, but it does not pinpoint it location. 

B.2. Superimposed Series of Flares 

Here we investigate the performance of the Maximum Peak Method applied to a light curve which consists of a series 
of superimposed flares. Flares 1 and 2 are examples. As an input time series we use the simulated light curve shown 
in Figure Al. 
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Figure B3 shows the result of the Maximum Peak Method applied to the series of superimposed flares. The method 
confirms that the time delay detected with Double Power Spectrum is consistent with the predictions and the method 
excludes significant ranges of possible time delays. 



Fig. B3.— Maximum Peak Method applied to the simulated light curve shown in Figure Al. 





